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ABSTRACT 

In this work we present cosmological TV-body simulations of the Local Universe with 
initial conditions constrained by the Two-Micron Redshift Survey (2MRS) within a 
cubic volume of 180 h~^ Mpc side- length centred at the Local Group. 

We use a self- consistent Bayesian based approach to explore the joint parameter 
space of primordial density fluctuations and peculiar velocity flelds, which are com- 
patible with the 2MRS galaxy distribution after cosmic evolution. This method (the 
KiGEN-code) includes the novel ALPT (Augmented Lagrangian Perturbation Theory) 
structure formation model which combines second order LPT (2LPT) on large scales 
with the spherical collapse model on small scales. Furthermore we describe coherent 
flows with 2LPT and include a dispersion term to model flngers-of-god (fogs) arising 
from virialised structures. These implementations are crucial to avoid artificial par- 
allel filamentary structures, which appear when using the structure formation model 
with the 2LPT approximation and the data with compressed fogs. We assume ACDM 
cosmology throughout our method. 

The recovered initial Gaussian fields are used to perform a set of 25 constrained 
simulations. Statistically this ensemble of simulations is in agreement with a reference 
set of 25 simulations based on randomly seeded Gaussian fiuctuations in terms of 
matter statistics, power-spectra and mass functions. 

Considering the entire volume of (180 h~^ Mpc)^ we obtain correlation coeflicients 
of about 98.3% for the cell-to-cell comparison between the simulated density flelds 
and the galaxy density fleld in log-space with Gaussian smoothing scales of rs = 
3.5 /i~^ Mpc (74% for rs = 1.4 /i~^ Mpc). The cross power-spectra show correlations 
with the galaxy distribution down to scales of 2.2 — 3.0 h~^ Mpc 

The simulations we present provide a fully nonlinear density and velocity fleld 
with a high level of correlation with the observed galaxy distribution at scales of a few 
Mpc. 
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1 INTRODUCTION 

The leading paradigm in cosmology is given by hierar- 
chical structure formation, in which the large-scale struc- 
ture we observe today formed through gravitationally 
driven subsequent mer^in^ of objects of increasing size (see 
IWhite k Reeslll978l : [FYv fc Peebledll978l l. 

Cosmological simulations have played a major role in 
confirming this theory, obtaining structures statistically sim- 
ilar to the observed ones (see e. g. Davis et al.l (Il985h: 
ISpringel et al.l (|2005r ) or for a recent review iKuhlen et al.l 
(|201^)). These studies use random seeded fluctuations per- 
mitting one only to compare theoretical models to obser- 
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vations based on statistical properties (see e. g. the recent 
work studying the likelihood of large structures IPark et al.l 
I2OI2I ). 

Such an approach relies on the cosmological principle 
implying that any sufficiently large observed region can be 
taken as a representative part of the Universe. Hence large 
surveyed volumes in which the Universe can be considered 
homogeneous are required to suppress the so-called "cosmic 
variance" . Large volume simulations are especially necessary 
to study the cosmic evolution of various signatures imprinted 
in the primordial density ffuct nations, such as baryon acous- 
tic oscillations (BAOs) (as inlKim etal .1120091 : IPrada et al.l 
I2OI2I : lAngulo et all I2OI2I : lAlimi et al.ll2012l ). However, nu- 
merical A'-body simulations become extremely expensive, 
when in addition to large volumes the surveying of faint 
objects is required (i. e. high resolutions). 
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cosmic variance 



A complementary approach to the 
problem consists of trying to reproduce observed structures 
in the Local Universe. This would allow for a direct compar- 
ison between theory and observations. In this case smaller 
volumes are sufficient, since the particular cosmic realisation 
of the Universe is taken into account. 

The concept of this approach is to reconstruct the ini- 
tial density fields compatible with the galaxy distribution 
given a particular structure formation model. Cosmologi- 
cal A/'-body simulations using the phase information cor- 
responding to these fields are expected to naturally re- 
produce the observed structures in the Local Universe. 
The quality of the phase information in the initial den- 
sity fields determines the accuracy with which structures 
are reproduced by a cosmological A/'-body simulation. These 
kind of simulations are usually called "constrained simu- 
lations" (CS). Various methods have been used to obtain 
Gaussian fie lds given a set of c onstraints (see th e sem i- 
nal work s by Bertschineer^ ('l987'l: 'Hoff man fc Riba'S (|l99lh : 
Ivan de W evgaert & Bertschinger (1996 ) and the more re - 
cent ones [ Kitaura &: EnfiUnl (l2008l): iKitaura et al.1 (|201Qh : 
IJasche fc Kitaural (|2010l l: lKitaura et al.1 (|20lJ VL 

Hitherto, a number of pioneering attempts have fol- 
lowed this approach showing the large number of difficul- 
ties one enco unters when tryin g to perform const rained sim- 
ulations fseelKolatt et^ll996 : Bistolas & Ho ffmanI Il998l : 
iMathis et aLll2QQ2l : iKlvpin et al...2003. : ■Lavaux.2oTol l . These 
challenges range from redshift-space distortions, shot-noise, 
galaxy bias, to finding the mapping between Eulerian and 
Lagrangian space. 

Constrained simulations of the local universe have been 
already extensively used as a numeri cal near-field cosmolog - 
ical laboratory (for an overview see iGottlober et al.l |2010|) . 
The predictions of Warm and Cold Dark Matter cosmolog- 
ical models have been compared with t he observed distri- 
bution of dwarfs in t he Local Universe (|Zavala et al.ll2009l : 
iTikhonov et al.l l2009l ) and the r eionisation history of the 
Local Group has been studied (|lliev et al.l 1 20 111 ) Zoomed 
high resolution simulations of the Local Group have been 
used to study the satellite s of the Local Group galaxies (see 
e. g. iKnebe et al . 2011; L beskind eralll201ll: iKnebe et al.1 



rrr 



I2OIII : Idi Cintio et al.ll201ll : |Di Cintio et al."201^ the for- 
mation of the Local Group ([For ero- Romero et al.l l201ll ) or 
recently the strippi ng of gas from dwarfs movi ng through 
the Local Web (see iBenitez-Llambav et al.ll2013l ). 

The methods applied so far, need however a deep revi- 
sion considering the recent theoretical advances and the ac- 
quisition of new data tracing the Local Universe. While the 
displacement field h as been neglected in the works using ve- 
locity data (see e. g. iKlvpin et al.ll2003l). which are unbiased 



real-space tracers (see iBertschinger et al.lll990l: Dekel et al.l 



I1990I : IZaroubi et al.' '1999": ICourtois et al.l 12012). displace- 
ments have been considered within linear Lagrangian per- 
turbation theory (|Zel'dovichl (|l970l ) app roximation) in most 
of the works usin g galaxy redshifts (jKolatt et al.l Il996l : 
iMathis et al.l I2OO2I ). Nevertheless, treating in the latter 
case coherent redshift-space distortions in a very approx- 
imate way and thus leading to squashed structures per- 
pendicul ar to the line of sight, as is apparent in the 
work of IMathis et al.l ((2002) • A more s ophist icated treat- 
ment has been employed in iLavauxl (|2010l ) within the 
Zeldovich approximation, however lacking an analysis of 



the real-space reconstruction. An improvement may be 
achieved within t his approximati on from velocity data as 
well, as shown in IPoumler et al.l (|2013l ). It has also been 
shown that higher order perturbation theory is significantly 
more precise than the Zeldovich approximation to estimate 
both the displacement and the velocity fields (see th e pi- 
oneering works of Gramann (1993); Monaco & Efstat hioul 
(1999); Narav anan et al.1 (|200ll ): iMohavaee et al.1 (|200M 
| 2006.) and more recen t works on t his su bject in 



Kitaura fc Angulol dJoTJ): iKitaura et al.1 (|2012l ): iNevrinckl 
(|2012l ): lKitaura k HelB* ('2012')). 

Previous works used inverse a pproaches, in which e i 
ther the smoothed dens ity fi e ld (iNusser k Dekell Il99! 
or the galaxies dPeeblesI 1 19891: Nusser &; Branching I2OO1 



Branchini et al.1 I2OO2I : iBrenier et al.1 l2003l : iMohavaee et aD 
2003b) are moved back in time according to the displace- 
ment field initially evaluated in Eulerian space, hereby ef- 
fectively assuming a one-to-one relation in the Eulerian- 
Lagrangian mapping. This mapping is degenerate after shell- 
crossing and a statistical approach becomes necessary. 

Alternative approaches have recently been proposed 
based on sampling the posterior distribution function of a 
Gaussian prior with a likelihood relating the initial fields 
with the observations through a particular structure for- 
mation model (see Kitaura 2013a; Jasche & Wandelt 20121: 
IWang et al.ll2013|) with its unprec edented application to ob- 



Kitaura et aD (|2012l ). 



ser vat ions in 

We present in this work for the first time constrained A/'- 
body simulations using such a statistical forward approach 
based on observations. In particular we aim at reproduc- 
ing the observed large-scale structure as it is traced by the 
distribution of galaxies in the Local Universe assuming hi- 
erarchical structure formation with ACDM cosmology and 
Gaussian seed fluctuations. 

We have introduced many improvements with respect 
to previous techniques by using the recently developed Aug- 
mented Lagrangian P erturbation Theory (A LPT) structure 
formation model (see IKitaura fc Hefil I2OI2I ) , and modelling 
redshift-space distortions taking into accou nt the tidal shea r 
tensor and including virialised motions fsee lKitaurall2013bl ). 

The quality of the constrained simulations which we 
present in this work permit us to reach a new level of pre- 
cision, in which we can effectively consider them fully non- 
linear reconstructions of the large-scale structure. The set 
of simulations we obtain is highly correlated on large scales 
(>98% with rs = 3.5/i~^Mpc), enabling us to study the 
influence of large scale modes on statistical processes (in a 
similar approach as Aragon Calvo 2013). The resemblance 
to the Local Universe makes it possible to directly compare 
the findings of the simulations with the observational data. 
The degree in which the observed distribution of matter is 
reproduced serves to test our cosmological models and verify 
our understanding of structure formation. A large variety of 
applications can be done based on these simulations with 
unprecedented accuracy, ranging from studies of the non- 
linear density and peculiar velocity fields, over galaxy bias 
analysis, to study merging histories and dependence on en- 
vironment of particular objects. 

The remainder of the paper is organised as follows. First 
we introduce the method, then we present the input data fol- 
lowed by our constrained simulations and finally we discuss 
the results and present the conclusions. 
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Figure 1. Methodology applied in this work. The left part of 
the flow-chart represents the KiGEN-code to flnd the large-scale 
primordial fluctuations. In an iterative self-consistent way the 
density fluctuations (^(^alpt) ^^^ sampled at Lagrangian coor- 
dinates {^alpt}- Structure formation is modelled within ALPT 
leading to the density fleld in Eulerian space S(xalpt)- Then 
redshift-space distortions (RSD) are added yielding (^(^jalpt)- A 
likelihood comparison (shaded circle) is done with the galaxy fleld 
{zg}- This process is iterated within a Gibbs- and Hamiltonian 
sampling scheme. The right part of the flow-chart shows the se- 
lection of samples according to their cross-correlation with the 
galaxy fleld and the subsequent generation of initial conditions 
(IC) for the A^-body runs from the white-noise flelds ?7alpt(9)- 
The results (6(x),6(z)) are analysed and compared with the ob- 
servations (shaded circle) to gain insights about our structure 
formation modelling and improve our reconstructions. 



2 METHOD 

Our method consists of exploring in a self- consistent 
Bayesian based approach the Gaussian primordial fluctua- 
tions which after nonlinear cosmic evolution are compatible 
with the galaxy distribution. Here, approximate and effi- 
cient gravity solvers are needed. We then scan the large set 
of reconstructions to select the ones which show the largest 
correlation with the data. The corresponding set of Gaus- 
sian fields is resampled on a finer mesh to generate the initial 
conditions, which we finally use to perform high resolution 
constrained A^-body simulations. These are employed to fur- 
ther confine the models used in the iterative scheme. Fig. [T] 
shows a flow- chart summarising the different steps in the 
method. 



2.1 Assumptions 

In our approach we assume that the primordial density fields 
are Gaussian distributed. Hence the two-point correlation 
function fully characterises the statistics of the initial con- 
ditions. _The_corresponding power-spectrum is derived with 
CAMB ( Lewis et a l. 2000) using the WMAP7+BAO+hO pa- 
rameters (|Komatsu et al.ll2011i ) . This cosmology is also used 



to evolve the density fields both in the reconstruction pro- 
cess using the approximate gravity solver (see ^2.2. ip and 
in the constrained simulations using an A^-body code. We 
assume that the density distribution in the Local Universe 
can be reproduced within the framework of ACDM cosmol- 
ogy and neglect for the time being baryonic physics. This 
will certainly play an important role at scales smaller than 
the Jeans length which is out of the scope of this work. Fur- 
thermore we assume periodic boundary conditions. This as- 
sumption is not critical, as the data considered here (see 33]) 
are increasingly sparse towards the boundaries (the selec- 



tion function drops to values of about 10 



10" 



While 



this will reduce the accuracy close to the boundaries, the 
inner regions are not affected considering that the correla- 
tion function of the density fluctuations quickly drops at 
distances larger than about 10 h~^ Mpc. Peculiar velocities 
on the other hand can be correlated up to distan ces of about 
100 h-^ Mpc (see e. g. lErdogdu fc Lahavll2009l ). Therefore, 
we will miss the external a ttractors t o the volume consid- 
ered in our work (see HoffmanI l200d : ICourtois et al.l I2OI2I , 
to add a posteriori the missing tidal component from ve- 
locity data) . The details of the different models used in the 
reconstruction process can be found below. 



2.2 Reconstruction of the large-scale primordial 
fluctuations 

We rely on the Ki GEN-code to rec over the large-scale primor- 
dial fluctuations (|Kitaurall2013al l . The method implemented 
in this code is based on sampling Gaussian fields constrained 
by Lagrangian test particles {q} (here Hamiltonian sam- 
pling is employed as developed in J asche fc Kitaural l20ld : 
iKitaura et al.|[201^ V which after cosmic evolution fsee l2.2.l"|) 
and once redshift-space distortions have been included (see 
I2.2.2P yields a distribution in Eulerian redshift-space {z} 
compatible with the observed galaxies {^^g}- The set of con- 
straints (matter tracers on a grid) in Lagrangian space {q} 
is obtained from the likelihood comparison between the par- 
ticles in Eulerian redshift-space with the galaxies, as the in- 
formation of each particle's trajectory is known. Both quan- 
tities, the density fluctuations and the set of particles in 
Lagrangian space, ar e iteratively obtained w ithin a Gibbs- 
sampling scheme fsee lKitaura fc EnBlinll 20081 , and references 
therein). Effectively we sample the posterior distribution 
function of Lagrangian density fields 6(q) given a set of 
tracers in redshift-space {^^g}, a structure formation model 
M and a peculiar velocity field v: P(S( q)\i zg} , M, v) ■ For 
more details on the method we refer to iKitaural (|2013a 



2. 2. 1 Structure formation model 

In order to scan the parameter space of primordial fluctua- 
tions we need an efficient and computationally fast gravity 
solver. Most approximate structure formation models are in- 
accurate when shell-crossing becomes dominant, i. e. when 
particles massively start crossing each others trajectories. 
Too severe shell-crossi ng can be prev ented with the spheri- 
cal collapse model f see lNevrinckll2012 ) . However, one needs 
to take into account the tidal field component to be accurate 
on large-scales (see e. g. Scoccimarro & Sheth 2002). There- 
fore, we use the ALPT approach (see Kitaura &: Hel3ll2012l V 
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This method is based on sphtting the displacement field into 
a long- and a short-range component. The long-range com- 
ponent is computed by second order LPT (2LPT). This ap- 
proximation contains a nonlocal and nonlinear tidal term. 
The short-range displacement field is given by a spherical 
collapse based solution. A Gaussian kernel with a smooth- 
ing radius of 4 h~^ Mpc is used to separate between both 
regimes. 



2.2.2 Redshift- space distortions treatment 

Coherent redshift-space distortions are modelled within 
2LPT. However, we consider two approaches to treat highly 
nonlinear redshift-space distortions. 

• CFOG: The first being the widely used one in 
which the fingers- of- god (fogs) are compres sed (see 
e. g. lErdogdu et ani2004l : Pfegmark et al]l2004al ). We use 
a friends- of -friends algorithm taking into account the ellip- 
soidal distribution along the l ine- of- sight of groups o f galax- 
ies due to virial motions (see Tegmark et al.ll2004 bV 

• SFOG: In t he second appro ach we apply a novel sam- 
pling sch eme fseelKitaurall2013bl). We follow the ansatz pro- 
posed in iKitaura fc EnBlinI ((20081) of splitting the velocity 
field into a coherent flow 1;^°^ and a dispersion term v^'. 
Vr — v%^ + Vr within the it erative reconstruction process 
(see also iKitaura et al.l I2OI2I ). The dispersion term mod- 
els the velocity distribution in closely virialised structures 
producing elongated structures along the line- of- sight. Here 
we use a density dependent dispersion and sample t he red- 
shift d i splacement fro m a Gaussian distribution fsee lKitaural 
l2013bl : lKitaural[2QQ7l ) . 

We will compare both approaches against each other in 
our numerical calculations section (Sec.[4j). 



2.2.3 Galaxy bias 

The galaxy bias which is effectively modelled within Ki- 

GEN includes a nonlinear, nonlocal and a scale-dependent 

component. In the likelihood comparison step haloes are 

constructed at the positions of galaxies within the ALPT 

approximation with a friends- of- friends scheme (for a re- 

I II 1 

view on the halo- model see lCoorav fc Sheth 2002). A simi- 

lar procedure was sug g ested for the 2LPT app roximation by 
IScoccimarro fc ShethI (|2002l ) and applied in iManera et al.l 
(|2012l ). Further, in the posterior sampling step of Gaus- 
sian fields a scale dependent bias is effectively modelled. 
The over-all shape of the power-spectrum is strongly de- 
termined by the prior correcting for a /c-dependent bias. 
We will show in ^4.31 that the characteristic features of the 
power-spectrum come from the particular matter distribu- 
tion in the considered volume. Sampling the full shape of 
the power-s pectrum within the reco n struction pro cess as 
suggested in IKitaura fc EnBhnl (|2008l ) : [jasche et al.1 (|2010 ): 
IKitaura et al.l ((20121 ) is out of the scope in this work. 



2.3 Scanning the posterior distribution function: 
selection of initial conditions 

Once the Markov chain has converged it starts producing 
Gaussian fields sampled from the posterior probability dis- 



tribution function (PDF). This enables us to compute vari- 
ous estimators like the mean of the posterior PDF y ielding 
a smoothed conservative field (see IKitaura et al.ll2010l . l2012l . 
and ^4.3p . We are, however, interested in samples which have 
the full power and at the same time can be considered as 
representative of the PDF. The mean properties of this sub- 
set should coincide with the corresponding mean of the full 
set having a smaller varianc e (see the case of peculiar mo- 
tions in .Kitaura et al.ll2012l ). Therefore we select the set of 
primordial fluctuations which show the highest correlation 
with the galaxy field after cosmic evolution within a par- 
ticular structure formation approximation. Here we choose 
a spectral analysis to estimate the degree of correlation be- 
tween the dark matter (DM) reconstructions (5dm and the 
galaxy overdensity field fc, i. e. the cross power-spectrum: 



XP(A:)[foM,fc] 



(|foM(fe)fc(fe)|) 



(1) 



where the ensemble brackets denote angular averaging and 
Pdm(A;) are Pg(A:) the power-spectra corresponding to the 
fields foM and fc, respectively. The matter overdensity field 
is defined by 



SbM : 



A^DM 



(2) 



where A^dmz denotes the number counts of simulated par- 
ticles in cell i gridded with a clouds-in-cell (GIG) scheme, 
and A^DM is the mean number counts computed dividing the 
total number of particles by the total number of cells in the 
considered volume. We define the galaxy overdensity by the 
maximum likelihood of a Poissonian distributi on to compen- 
sate for selection function effects /s (see e. g. IKitaura et al.l 
,2009,): 






1. 



(3) 



for each cell i, with the mean expected number counts being 
given by: A'g = J2i^Gi/fsi- We impose a threshold XPth 
at a given scale to filter out the tails of the posterior PDF 
and focus on the denser probability regions. Approaching 
the grid-resolution at high k the cross power-spectra start 
to fiuctuate. Therefore we choose a value of kth far below 
the Nyquist frequency (kth — 1, see ^4.2p : 

^^kLPT,G(A: = /Cth) = XP\k = /Cth)[^ALPT,fc] > APth , 

(4) 
for sample / in the Markov chain. 

Gross power-spectra which satisfy this condition usu- 
ally show higher correlations over the entire /c-space region 
(see M.2p . To avoid selecting correlated samples, which are 
essentially identical, we discard first and second neighbours 
taking the largest correlated one among them: 

^^alpt,g(^ = ^th) > max(AP5^LPT,G^^(^ = ^th)) . (5) 



2.4 Generation of high resolution initial 
conditions 

From the set of primordial fiuctuations we construct the 
corresponding white-noise fields: 



ri\k) 



5\{k) 



(6) 
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for sample /, with P{k) being the theoretical linear power- 
spectrum and Sq standing for the overdensity field at the 
initial conditions in Lagrangian space q. We note that these 
fields are essentially uncorrelated from cell-to-cell, except for 
some residual deviations from a flat power-spectrum char- 
acteristic of the particular cosmic realisation. We checked 
that dividing Sq{k) by the square root of the actual power- 
spectrum measured from the density field itself, i. e. with 
perfectly flat white-noise fields, leads to significantly less 
correlated density fields with the galaxy field after simu- 
lating structure formation. The construction of the white- 
noise field permits us to increase the resolution augment- 
ing the smaller scales on a finer mesh randomly sampling 
from a Gaussian with unity variance. The resulting fields can 
then be multiplied in Fourier-space with the linear power- 
spectrum including modes up to the resampled scales. Here 
we assume the same cosmology as in the reconstruction al- 
gorithm to maximise the precision of the constrained sim- 
ulations. From the resulting density fluctuations one needs 
to compute the displacement and velocities for the number 
of particles required in the A/'-body simulation. These steps 
are done within the initial conditions (ICs) generation code 
GINNUNGAGAfIj, which is fully parallelised and has been ex- 
tensively tested within the Zeldovich approximation. 



2.5 Constrained A/^-body simulations 

The initial conditions computed as described in the previous 
sub-section are used to simulate structure formation with a 
iV-body code. To avoid transients which could originate from 
the Zeldovich approximation in the generation of the ICs, we 
start the A/'-body simulations at a relatively high redshift of 
z = 100. We use GADGET-3, an improved version of the pub- 
lic ly available co smological code gadget- 2 (last described 
in ISpringell l2005l ) . The code is parallelised for distributed 
memory machines and uses a hierarchical tree algorithm for 
the calculation of gravity. We will restrict ourselves in this 
work to dark matter only A/'-body simulations and leave a 
detailed description of the star-forming, cooling and feed- 
back processes for future work. The advantage of this step 
is that we can accurately model the nonlinear regime and 
get a realistic distribution of haloes. Haloes are defined in 
this context as spherical regions with a density higher than 
p > 200 X pcrit, where pcrit is the critica l density. To identify 
these regions we use the AHF-codqj ( Knollmann fc Knebd 
I2OO9I I . The halo distribution can be used to further test the 
reconstruction algorithm and tune the models as we will 
show below (see ^221 and 34T]). 



3 INPUT DATA 

Our input data in this work is the 2MASS redshift sur- 
vey (2MRS), Ks = 11.75 catalogu e, with a median redshift 
of z ^ 0.028 (jHuchra et al.ll2012l ). This catalogue provides 
an almost full sky observation of galaxy redshifts with a 
sky coverage of 91% and uniform completeness of 97.6% 
only limited by the presence of the foreground stars near 



^ publicly available at https://code.google.eom/p/ginnungagap/ 
^ publicly available |http://popia. ft. uam.es/AHFJ 



the Galactic plane (the Zone of Avoidance (ZoA), where 
\b\ < 5° and \b\ < 8° near the Galactic centre). The same 
version of the augmented 2MRS catalogue was used as in 
I Kit aura et al.l (|2012l l, which was filled in the ZoA with 
random galaxies generated from the cor responding long i- 
tude/distance bins in the adjacent strips ( Yahil et al .11 19911 ). 
This method is robust for the width of the 2 MRS mask and 
has been thoroughly tested (for details see lErdogdu et al.l 
I2OO6I . and references therein). 

The ZoA could be sampled with a Poissonian likeli- 
hood describing the counts-in-cells of the galaxy distribu- 
tion wlnch_wouldb£2iipited to a coarse grid resolution (see 
e.g. iKitaura et al.ll201Ql ). To incorporate a mask treatment 
within a particle based reconstruction method like kigen, 
we would need to sample mock galaxies (or haloes) accord- 
ing to our structure formation model following schemes like 
the one proposed in Scoccimarro fc Sheth (2002) , which is 
out of scope in this work. 

We computed the selection function /s directly based 
on the rad ial distribution galaxie s ignoring the Kaiser-rocket 
effect f see JBranchini et al.ll2012i ). 



4 NUMERICAL CALCULATIONS 

In this study we focus on the Local Universe as comprised 
by a box of 180 h~^ Mpc on a side with the observer located 
in the centre including 31017 2MRS galaxies (i. e. a number 
density of 5.3 10"^ (/i"^ Mpc)"^). 

First we apply the KiGEN-code to find the large-scale 
primordial fluctuations. The resolution of the reconstruc- 
tion is about 1.4/i~^ Mpc on a mesh of 128^ cells. We follow 
iKitaura et al.l (|2012[ ) for the numerical setting of the KIGEN- 
code. 

The number of constraints (=total number of cells 
^ 2 10^) used in the Hamiltonian sampling process (see 
Sec J2.2p is achieved by taking into account all volume el- 
ements including empty cells. Therefore we account for ob- 
served regions with no detections. This amount of informa- 
tion excels that of previous works which relied exclusively on 
the locations of density peaks (see Sec. [T)) by about 2-3 or- 
ders of magnitudes. Even if one considers that peak positions 
are less correlated constraints than gridded volume elements, 
the latter obviously carries more information. The phase in- 
formation is used without any additional smoothing to the 
griding (see IWang et al.l l2013l , for smoothed constraints) . 
This permits us to recover not only the density peaks, but 
also filamentary structures as traced by the galaxy distribu- 
tion (see ^42]) . 

The power-spectra of the subsequent primordial density 
fields in the Markov chain are unbiased with the theoretical 
linear power-spectrum on all scales as soon as the kigen- 
code has converged. For safety we consider samples after 
2000 iterations for our study well after the previous criterion 
has been fulfilled. The Markov chain continues running to 
obtain a minimum of 1000 samples. 

Our reconstructions based on 2LPT and the ALPT ap- 
proximation using data with collapsed fogs (2LPT+CF0G 
and ALPT+CFOG, respectively) lead to constrained simu- 
lations with pronounced artifacts. These are caused by shell- 
crossing (see ^4.ip . As a consequence, the mass function of 
the case 2LPT-CF0G clearly deviated from the expected 
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Figure 2. Sketch of overlap- 
ping fogs (light grey) close to 
the observer at the bottom. 
A FOF based fog compression 
might accidentally detect the 
darker area and combine the 
fogs of separated objects (cen- 
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one object closer to the ob- 
server (light cross). See the 
case of the Virgo cand. (Fig.O 
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Figure 3. Artifacts due to shell-crossing and fog compres- 
sion (CFOG). Logarithmic DM overdensity in real-space in 
a slice of 22.5 h~^ Mpc thickness in the supergalactic XY 
plane. The panels correspond to a constrained simulation with 
(top:) 2LPT+CFOG, (middle:) ALPT+CFOG and (bottom:) 
ALPT+SFOG. The yellow ellipses show the negative fogs. The 
green highlighted area shows the accordion-effect, a preferential 
filament orientation perpendicular to the line-of-sight due to shell- 
crossing. The position of the Virgo cluster candidate (cyan circle) 
is shifted (2-3 h~^Mpc) towards the observer (red cross) in the 
upper two panels w.r.t. the case shown in the lower panel. 
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Figure 4. Left panel: PDF of finding the Virgo cluster candi- 
date in the 25 constrained simulations (ALPT+SFOG) at a cer- 
tain distance with respect to its mean position {x). Right panel: 
corresponding mass PDF from the 25 simulations. 



one based on random simulations. Therefore, we focus in 
this work on reconstructions using the ALPT structure for- 
mation model with the SFOG (ALPT+SFOG) prescription. 
We apply the full machinery sketched in Fig. [1] to perform 
constrained simulations. The simulations are conducted with 
384^ particles in a simulation box of 180 h~^ Mpc per side 
which corresponds to a mass of 7.8 x 10^ h~^ Mpc of each 
mass element in the chosen cosmology. To this end we used 
a comoving gravitational softening of 15 h~^ kpc. 

Let us assess the quality of the constrained A/'-body 
simulations. We will first analyse our redshift-space distor- 
tions (RSD) treatment followed by a comparison with the 
observed galaxy field. Finally, we will compare the statistical 
properties of the CSs with a set of random simulations. 



4.1 Redshift-space distortions treatment 

To investigate the accuracy of our SFOG treatment we 
compare it w it h the classical CFOG scheme following 
[Tegmark et al.l (|2004al ). The latter approach shows sig- 
nificantly larger correlations than the SFOG one. While 
the CFOG samples reach cross power-spectrum values of 
XP{k = 1)^0.7, the SFOG ones remain at XP{k = 
1)^0.5. The artificial collapse of fogs in the CFOG approach 
leads to relatively more clustered regions. This partially 
compensates for shell-crossing in the final density fields, and 
the corresponding lack of power of the approximate struc- 
ture formation model. This effect can induce higher cor- 
relations. However, a deeper investigation, shows that the 
CFOG approach induces artifacts. 

The first one comes from unfolding redshift-space to re- 
cover the real-space matter distribution. This happens in 
KIGEN in a forward approach, by adding redshift-space dis- 
tortions to th e simulated par ticles in the likelihood compar- 
ison step fsee lKitaurall2013al ) . A correct modelling of RSD is 
therefore necessary to obtain accurate reconstructions. One 
of the problems of using an approximate structure forma- 
tion model based on perturbation theory is that the highly 
nonlinear virialised motions are not included. As shown in 
iHamiltonl (|l998h coherent redshift-space distortions can al- 
ready produce fogs. Galaxies coherently approaching each 
other in a potential well, would produce shell-crossing in 
redshift-space due to their peculiar motions, leading to elon- 
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Figure 5. Cross power-spectra XP(k) between the DM reconstructions and the galaxy field. Left panel: XP(k) of a set of 1200 
realisations corresponding to the ALPT+SFOG reconstructions. The shaded area denotes the standard deviation of the complete set. 
The coloured lines show a subset of 25 realisations that have been picked according to criterion O Middle panel: Cross power-spectra 
comparing the 2MRS galaxy overdensity with the ALPT overdensities of the 25 reconstructions (solid) and for their corresponding 25 
constrained simulations (dotted). Right panel: Ensemble mean of the cross power-spectra between the galaxy field and the estimated 
dark matter field corresponding the ALPT reconstructions in redshift-space (solid green), the constrained simulations without RSD 
(dashed purple), with DM velocity based RSD (dotted blue) and with RSD based on Halo velocity fields (dash-dotted red). 



gated structures. Since in the CFOG case kigen tries to 
match the compressed data, negative fogs appear in real- 
space. This effect is boosted with artificial shell-crossing, as 
can be seen when using the 2LPT approximation (see fogs 
in real-space in the upper two panels of Fig. O where we 
zoomed into an area of 80 h~^ Mpc). Treating only coherent 
flows also fans out parallel filamentary structures like an 
accordion, as is apparent in the green highlighted areas in 
Fig. O One can see in the extension of the parallel filamen- 
tary structures that shell-crossing is more severe with 2LPT 
than ALPT (compare green areas). Another important as- 
pect is the uncertainty in the location of the mass distribu- 
tion in high density regions. While the CFOG scheme ignores 
this point yielding a single solution, the SFOG sampling 
approach propagates the uncertainty in each Markov chain 
sample. A fogs compressed based reconstruction is there- 
fore underestimating the uncertainty in high density regions. 
Moreover, the CFOG compression produces systematic arti- 
ficial displacements of clusters located close to the observer. 
This situation is depicted in Fig. O In regions where var- 
ious fogs overlap it is difficult to find the centres of mass 
with friends-of- friends algorithm. We found that clusters in 
such fogs are systematically located closer to the observer. 
Galaxies residing in the distant tails of fogs, which are more 
separated, may not be considered friends and hence a sin- 
gle shorter fog is detected with its centre closer to the tails 
where the fogs converge, i. e. closer to the observer. This 
is the case of the local super cluster which includes Virgo. 
Fig. [3] shows the location of our Virgo candidate (massive 
cluster within the cyan circle) with respect to the observer 
(red cross) with both schemes. We find that the position of 
the Virgo candidate is estimated to be about 2.5/i~^Mpc 
closer to the observer with the CFOG scheme than with the 
SFOG one (distances of 6.8 h~^Mpc and 9.3 h~^Mpc respec- 
tively for the ensemble means of 25 simulations each) . 



The robustness of the RSD treatment can be seen in 
the uncertainty of the position of particular objects. The 
spread in distance of the Virgo candidate across all 25 
ALPT+SFOG constrained simulations is very small (1 a 
< 1 h~^ Mpc), as shown in Fig. (4] Also the mass of the cor- 
responding object is strongly constrained as shown in the 
right panel. This indicates that a constrained environment 
leads to a constrained final mass of halos. We plan to inves- 
tigate this further and plan to study the mass evolution of 
particular objects and the influence of the environment. 

It is possible to improve the SFOG scheme in various 
ways. The description of coherent velocities should be re- 
vised. Furthermore, a radial dependent dispersion velocity 
needs to be implemented to account for observational se- 
lection effects, as we plan to do in future work. Neverthe- 
less, the SFOG approach as presented in this work leads to 
constrained simulations, which solve many of the problems 
present with the common CFOG scheme to great extent. For 
this reason we will restrict from now on our study to the set 
of 25 ALPT+SFOG based constrained simulations. 



4.2 Comparison with the observations 

The cross power-spectra between the ALPT reconstructions 
with the SFOG treatment and the galaxy field are shown 
on the left panel in Fig. \5\ They present very high values 
up to /c ~ l/iMpc~^ and are still correlated at 2.2 < k < 
3 /iMpc~^ (1 a area of the zero crossing of the cross power- 
spectra). We select the 25 most correlated samples repre- 
sented by the different coloured lines as shown in Fig. [5l 
Those samples, which are higher correlated with the galaxy 
field at a given /c, show in most of the cases higher correla- 
tions in the entire /c-space. Such a rank-ordered behaviour 
is broken in some cases when approaching the Nyquist fre- 
quency at k>2hMpc~^. This demonstrates the robustness 
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of our selection approach focusing at /c = l/iMpc~^ (see 
^2.3p . The required threshold at /c = 1 to select the 25 sam- 
ples was: XPth = 0.43. The middle panel of Fig. \5\ shows 
the cross power-spectra corresponding to the A/'-body simu- 
lations, being systematically less correlated as compared to 
the ALPT ones. The reason for this is that the initial con- 
ditions have essentially been computed within the ALPT 
approximation in addition to a particular peculiar velocity 
model (2LPT with a dispersion term) to simulate RSD. The 
iV-body solutions (modelling gravity more accurately) do 
not fully agree with these models and introduce a small dis- 
placement. It would thus be interesting to include in the 
future full gravity solvers within the reconstruction process. 
In the right panel we demonstrate that our RSD treatment 
is over-all correct yielding higher correlations when trans- 
forming the simulation from real-space (violet curve) into 
redshift-space (cyan and red curves). Here we consider two 
ways of transforming the data, using the peculiar velocities 
of the dark matter particles in one case (cyan) and the ones 
of the halo in the other (red). The right panel in Fig. [5] 
demonstrates that our RSD modelling is over-all correct as 
the cross-power spectra are clearly higher when including 
redshift-space distortions. The lower correlation when us- 
ing the halo velocity fields is due to the additional stochas- 
tic process of collapsing objects along the line-of-site. The 
larger dispersion when taking the DM velocities enhances 
the probability that a number of simulated particles match 
the positions of observed galaxies. Nevertheless, matter is 
distributed in a more realistic manner using the halo ve- 
locity fields as we will discuss below. We would like to em- 
phasise that we take the galaxy field according to Eq. [3] 
including redshift-space distortions (also fogs) in our cross- 
correlation analysis (hence performin g a direct compar ison 
with the galaxy density field, where Wang et al.l l2013l , for 
example, cross-correlate with a processed field). 

To further assess the quality of the constrained simula- 
tions we make a cell-to-cell comparison between the simu- 
lated and the observed density fields in configuration space. 
The correlation between two fields ^i and 62 can be quanti- 
fied through the Pearson correlation coefficient r 



r(6i,62 



E^(^l^-^l)(^2^-^2) 



(E.(^H-^i)) (E,(^2,-^2))^ 



(7) 



Let us focus here on the logarithmic overdensity ln(l + 5) for 
various reasons: first this representatio n yields an estimate of 
the linear density field c omponent f see iNevrinck et al]|2009l : 
iKitaura &: Anguldl2012l l: second it suppresses shot-noise in 
overdense areaqj. Fig. [6] shows the average cell-to-cell cor- 
relation between the A/'-body simulations and the galaxy 
field, and respectively ALPT, in different grid resolutions as 
compared to the ALPT reconstructions. We find that the 
uncertainty is much larger in the underdense regions, as can 
be seen in the increasing scatter towards negative values. 



^ Let us write the galaxy overdensity field as the sum of the ex- 
pected galaxy density field and the shot-noise term: ^g = ^g^^ +^ 
(see Kitaura et al. 2010). The residual shot-noise term is accord- 
ingly given in the logarithmic representation by: ln(l + (5g) — 
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Figure 6. Two dimensional histogram of a cell-to-cell com- 
parison in logarithmic-space including RSD between the con- 
strained simulations and the (top:) 2MRS galaxy field; (bot- 
tom:) ALPT-hSFOG reconstructions. A CIC algorithm was 
used with (left:) 64^ cells and Gaussian smoothing with rs = 
3.5/1-1 Mpc; (right:) 128^ cells and rs = 1.4 /i^^ Mpc (cell 
width). Darker grey-scale indicates a larger number of cells. 



The close alignment with the 45° slope in the upper pan- 
els demonstrates that the galaxy bias is very close to unity, 
especially on scales larger than rs = 3.5 h~^ Mpc. The de- 
viations from the 45° slope shown in the lower panels are 
due to the perturbative approach of the ALPT approxima- 
tion, which does not fully cap ture the nonlinear be haviour 
of gravitational clustering f see iKit aura &; Hedl2012l ). 

The corresponding cross-correlation coefficients are very 
high, being r = 98.3 ±0.1 % on smoothing scales of 
rs = 3.5 h~^ Mpc for both the ALPT reconstructions and 
the constrained simulations. For rs = 1.4 h~^ Mpc we get 
r = 78.6 ± 0.3 % and 73.8 ± 0.1 %, respectively. It is 
remarkable that we still obtain r = 97.5 ±0.1 % with 
rs = 3.5 h~^ Mpc for the cross-correlation between the 
galaxy field (in redshift-space) and the constrained simu- 
lations in real-space. This tells us that small changes in r 
of ^ 1% with rs = 3.5 h~^ Mpc can imply improvements in 
the precision of the reconstruction of the order of ^-^ 5h~^ 
Mpc (the order of coherent RSD). Such an estimatio n is im- 
porta nt when comparing with previous works (see iLavauxl 
l20ld where a value of r = 97% is obtained taking a smaller 
volume made by a sphere with a radius of 60 h~^ Mpc and 
a volume limited 2MRS galaxy field using the same rs). 

Fig. [71 visualises the accuracy of the constrained simu- 
lation by plotting on top of the dark matter density field 
the distribution of galaxies. It is remarkable how accurately 
the red dots, each one standing for a galaxy, match the sim- 
ulated structures in redshift-space. The constrained simu- 
lations are even able to connect filaments traced by a few 



© 0000 RAS, MNRAS 000, 000-000 



Structure Formation of the Local Universe 9 



I 






I 






I 



o 




50 



Q 
^ 



-50 



50 



Q 

^ 



-50 



50 



en 
O 

^ 



-50 



50 -50 

SG X [/i-i Mpc] SG Z [/z-i Mpc] 

Figure 7. Density fields in supergalactic coordinates. Upper panels: redshift-space. Logarithmic DM overdensity of one of our 
constrained simulations (based on ALPT-hSFOG) in a slice of 22.5 h~^Mpc thickness in supergalactic coordinates XY around Z=0 (left) 
and ZY around X=0 (right). Redshift-space distortions have been computed according to the DM particle velocities. The overdensity 
maps have been blended with the distribution of galaxies (red dots). Middle panels: redshift-space. Same as upper panels, but taking 
the halo velocity field. Lower panels: real-space. Same as upper panels without galaxies. The densities have been calculated with a 
SPH-kernel fixed to 32 neighbours. 
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Figure 8. Signal and signal-to-noise maps in the supergalactic XY plane. Left panels: Ensemble-averaged logarithm of the overdensity 
fields {\og(p/p)) of ALPT+SFOG reconstructions (top) and constrained simulations (bottom) showing the middle slice of a 128^ CIC- 
grid. Right panels: Ensemble-averages divided by the standard deviation of the ensemble of reconstructions (top) and simulations 
(bottom) see Eq. [8] 



galaxies. The upper panels show the maps using the DM ve- 
locities, while the middle ones are based on the halo velocity 
field. For definiteness we move the simulation particles ac- 
cording to the velocity of their host halo if available. For non 
virialised particles we use the halo velocity field on a grid 
of 128^. We can clearly see that in the upper panels fogs 
are overestimated. Structure formation forms bound objects, 
therefore the correct way of modelling RSD consists of tak- 
ing the halo velocity field to transform the simulation from 
real-space to redshift-space. To accurately model the RSD 
one needs a high enough resolution to resolve the masses of 
the observed objects and consider selection effects. The ve- 
locity dispersion of galaxies decreases with distance to the 
observer, since magnitude cuts tend to select more massive 
objects at increasing distances. We plan to perform a study 
of high resolution A/'-body simulations to study this in de- 
tail. Our current simulations yield already realistic fogs, as 



can be seen from the middle panels in Fig. [71 The lower 
panels in that figure show a natural real-space matter distri- 
bution with more extended structures along the line-of-site 
as opposed to the squashed structures in redshift-space (see 
Coma, Great Attractor and Perseus-Pisces region). Here we 
also find that the Local Void is less empty than it appears 
in redshift-space. 

4.3 Statistical analysis 

Here we want to investigate the robustness of the con- 
strained simulations by looking at their statistical proper- 
ties. Fig. \8\ shows in the left panels the ensemble- averaged 
densities in logarithmic space (log(^ + 1)) = (log {p/~p)) cor- 
responding to the set of ALPT reconstructions and simula- 
tions. This demonstrates the consistency between the ALPT 
and the A/'-body simulations on scales larger than a few Mpc 
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Figure 9. Matter statistics. Left panel: Matter PDF of the 25 
constrained simulations (same colour-code as Fig. [5]) at initial 
redshift z = 100. Right panel: same as left panel but at final 
time z = 0. The light and dark grey shaded contours indicate 1 
and 2 sigma regions according to the 25 random simulations. 



confirming our cell-to- cell comparison. However the A/'-body 
simulations show as expected, much more small-scale power. 
Looking at the signal-to-noise ratio: 



S_ 

N 



(log(l + foM)) 



^((log (1 + foM))' - (log (1 + foM))'> 



(8) 



where the ensemble mean is computed by volume averag- 
ing, we find that the uncertainty is smaller (darker colours) 
for the ALPT case. The small-scale power in the A/'-body 
case can be especially seen in the larger uncertainties in low 
density regions. 

To accurately assess the statistical quality of the con- 
strained simulations we construct a reference sample of 25 
randomly seeded A^-body simulations with the same setup. 

We checked that the white-noise fields used for the ini- 
tial conditions have numerically vanishing statistical mo- 
ments: mean, skewness and kurtosis. Fig. [9] shows the mat- 
ter statistics for the constrained simulations at redshifts 100 
and 0. We find that the curves lie within the variance of the 
reference random simulations. The mass function (Fig. llOp 
shows a slight overproduction of massive haloes for the mass 
range M200 > 10^^ /i^^M©, being, however, compatible with 
the reference set within the error bars. It is unclear, whether 
the particular realisation of the observed volume could be 
responsible for this issue. 

Finally, we analyse the power-spectra to verify the 
method and the quality of the constrained simulations. This 
is shown in Fig. 1111 Here one can see that the specific fea- 
tures of the galaxy power-spectrum, characteristic of the 
particular matter distribution, are preserved both in the 
ALPT reconstructions and in the constrained simulations. 
We can also see that the CSs have more power than the 
ALPT reconstructions, correcting for the bias caused by the 
perturbative structure formation approach. 

The black solid line shows the excess of power for the 
galaxy field due to the galaxy bias and shot-noise. The solid 
lines show the strongly biased nature of the ALPT+SFOG 
model. The constrained simulations however, are close to 
unbiased with respect to the random simulations and show 
the characteristic features of the particular realisation of the 
Local Universe. Fig. [TT] demonstrates that we are reasonably 
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modeling, in a scale dependent way, the biased nature of the 
galaxy field and the nonlinear gravitational clustering up to 
high ks {k^2hMpc~^). 



5 CONCLUSIONS 

We have presented in this work constrained A/'-body sim- 
ulations of the Local Universe based on the 2MRS galaxy 
redshift catalogue. This is the first time that a statistical 
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self- consistent forward approach is applied to observations to 
obtain the fully nonlinear phase-space distribution. Our ap- 
proach is based on sampling the posterior distribution func- 
tion of a Gaussian prior with a likelihood relating the initial 
fields with the observations through a particular structure 
formation model. In this way we can exploit the high num- 
ber density of observed galaxies in the Local Universe and 
use between 2 and 3 orders of magnitude more constraintsfor 
the initial fields than previous works. Special care has been 
taken to include recent improvements in modelling structure 
formation and redshift-space distortions within the recon- 
struction method of the primordial fluctuations. In partic- 
ular, Augmented Lagrangian Perturbation Theory (ALPT) 
was used and redshift-space distortions were modelled, hence 
taking into account the tidal shear tensor and including the 
modelling of virialised motions. 

Our studies have shown, that this effort is necessary to 
deal with such a large number of constraints, while resolv- 
ing the density field down to scales of 2-3 /i~^Mpc. Other- 
wise, strong artifacts are introduced in the constrained sim- 
ulations due to transients which originated through shell- 
crossing in the iterative reconstruction process. We verified 
that this is the case, when using standard second order 
Lagrangian perturbation theory or when using the galaxy 
distribution where the fingers-of-god have been collapsed. 
Parallel filamentary structures perpendicular to the line-of- 
sight and negative fingers-of-god in real-space appear in such 
cases. We also found that the position of clusters which 
are distorted as fingers-of-god close to the observer (such 
as Virgo) tend to be systematically closer in the collapsing 
schemes. We have managed to strongly suppress all these 
effects with our improved treatments. Nevertheless, we con- 
sider that further investigation should be carried out im- 
proving the models with respect to the statistical properties 
of the constrained simulations. One focus should be on the 
mass functions. We found, that those show a slight overpro- 
duction of massive haloes, which may or may not be induced 
by an artifact. This is not clear, as the deviation is small and 
can be considered compatible with the random simulations. 
Especially so since one must also consider that the Local 
Universe is a particular realisation, which could have very 
specific statistical characteristics. Nevertheless the compat- 
ibility between the constrained simulations and the random 
ones in the matter statistics at starting and final redshifts 
and in the power-spectra, demonstrates the quality of the 
constrained simulations. Certainly, we should improve our 
selection function treatment, including a real- space estima- 
tion to treat the Kaiser-rocket effect. We have considered 
that this is negligible given the small volume we are study- 
ing. In any case, such an improvement would also permit 
us to tackle larger volumes for which the selection function 
drops more dramatically. In those studies one could also 
investigate the peculiar motions on larger scales, including 
more distant attractors, since now we are limited by the size 
of the box (180 /i~^ Mpc) and the periodic boundary condi- 
tions. 

This work serves as a consistency check for our under- 
standing of structure formation. In this sense the high level 
of resemblance between our models and the actual galaxy 
distribution corroborate the initial assumptions and there- 
fore do not indicate any deviations from the standard cos- 
mological model. It would be interesting to extend this kind 



of study varying the assumptions to further test our models 
and verify whether these are the most compatible with ob- 
servations. The level of precision that we achieve already in 
this work should help to get new insights into fields of study 
that require accurate density estimates of the Local Universe 
such as cosmic rays, the warm hot inter-galactic medium or 
Dark Matter annihilation and decay signals. The full non- 
linear velocity field determined in this work can contribute 
to our understanding of the cosmic flows in the Local Uni- 
verse. Furthermore the simulations can also provide a con- 
sistent formation history of particular clusters. We plan to 
exploit this by studying objects within their correct environ- 
ment, explore their formation histories and compare them 
with observations. 
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